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Abstract A numerical model is used to simulate rheome¬ 
ter experiments at constant normal stress on dense sus¬ 
pensions of spheres. The complete model includes sphere- 
sphere contacts using a soft contact approach, short 
range hydrodynamic interactions defined by frame-invariant 
expressions of forces and torques in the lubrication ap¬ 
proximation, and drag forces resulting from the porome- 
chanical coupling computed with the DEM-PFV tech¬ 
nique. Series of simulations in which some of the cou¬ 
pling terms are neglected highlight the role of the porome- 
chanical coupling in the transient regimes. They also re¬ 
veal that the shear component of the lubrication forces, 
though frequently neglected in the literature, has a dom¬ 
inant effect in the volume changes. On the other hand, 
the effects of lubrication torques are much less signifi¬ 
cant. 

The bulk shear stress is decomposed into contact 
stress and hydrodynamic stress terms whose depen¬ 
dency on a dimensionless shear rate - the so called vis¬ 
cous number I v - are examined. Both contributions are 
increasing functions of / v , contacts contribution domi¬ 
nates at low viscous number (I v <0.15) whereas lubri¬ 
cation contributions are dominant for I v > 0.15, consis¬ 
tently with a phenomenological law infered by other au¬ 
thors. Statistics of microstructural variables highlight 
a complex interplay between solid contacts and hydro- 
dynamic interactions. In contrast with a popular idea, 
the results suggest that lubrication may not necessarily 
reduce the contribution of contact forces to the bulk 
shear stress. The proposed model is general and applies 
directly to sheared immersed granular media in which 
pore pressure feedback plays a key role (triggering of 
avalanches, liquefaction). 

Keywords granular suspension • rheology • lu¬ 
brication • shear flow • discrete element method • 
hydromechanical coupling 
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1 Introduction 

Dense suspensions of particles immersed in a viscous 
fluid are ubiquitous in natural phenomenon, such as 
sediment transport or debris flows, and in numerous in¬ 
dustrial applications such as form filling with fluid, con¬ 
crete in civil engineering or slurry transport in petroleum 
industries. The understanding of dense suspension rhe¬ 
ology has lead to an important research effort over the 
past decades pifTMT]. The complexity of this problem 
arises from its two-phase nature involving a fluid phase 
(continuous) and a particulate phase (discrete) for which 
particle-particle interactions and fluid-particle interac¬ 
tions contribute to the behavior of the system in the 
dense limit. 

Classical rheometer experiments impose simple shear 
of the suspensions at constant volume (type I). In such 
case, the interpretation of the shear stress in terms of 
effective viscosity r\ e suggests that rj e diverges at high 
solid fraction (</> « 0.6 for spheres) [HI). Rheometer ex¬ 
periments at constant normal stress P have been per¬ 
formed only recently [5] (type II). Under such condi¬ 
tions the volume of the suspension is free to change 
as a response to the imposed shearing. On this basis, 
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a description of dense suspensions has been proposed, 
which unifies classical suspension rheology, described 
in terms of a shear and normal effective viscosity which 
depend on the solid fraction, and the dense granular 
flow rheology, described in terms of shear to normal 
stress ratio (/i) and solid fraction (</>). For this purpose 
a dimensionless shear rate was introduced, the so called 
viscous number (I v ), controlling both frictional and vis¬ 
cous contributions to the shear stress. An advantage of 
such a visco-plastic vision is that the behaviour of the 
suspension in the very dense limit is not associated with 
a divergence of the viscosity. Instead, the shear to nor¬ 
mal stress ratio reaches a constant value corresponding 
to a Coulomb-type bulk friction. 

From an analytical point of view, the rheology of 
suspensions has been studied since the beginning of the 
20st century. Einstein in 1905 m derived the effective 
viscosity of a dilute suspension based on long range hy¬ 
drodynamic interactions. Frankel and Acrivos m pro¬ 
posed another derivation in the limit of dense packings. 
The later was based on the so-called lubrication ap¬ 
proximation of the hydrodynamic interactions between 
nearly touching spheres. The lubrication terms are sin¬ 
gular, they appear as pair interactions between particles 
and they diverge at the approach of contact. 

The first discrete numerical simulations of parti¬ 
cle suspension has been proposed in the framework of 
Stokesian Dynamics (SD) [4,6] using resistance and 
mobility matrices [18,. This technique is able to quanti¬ 
tatively reproduce the divergence of effective viscosity 
for solid fraction approaching close packing. In the gen¬ 
eral framework of SD, the hydrodynamic forces include 
both long range and short range interactions which are 
defined independently, the later are the diverging terms 
as found in the lubrication approximations. In addition 
the forces depend on fluid velocities and particle veloci- 
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ties as independent variables. Practical implementation 
of this general framework in computer codes for many 
particles is still a great challenge, however. Commonly, 
the fluid velocity unknowns are eliminated by assum¬ 
ing that the fluid comoves with the particles at large 
scales [H[6l fCT[T9l[2Ql[TQl [3] . This assumption also entails 
that the pair drag forces must be frame invariant [3], 
which exclude many components of the full resistance 
matrices. In fact, all the long range interactions must 
be removed, and only the lubrication terms are left. 

The assumption of co-movement is acceptable for 
shear flow at constant volume (type I), but it is oth¬ 
erwise a severe limitation of the numerical technique. 
One consequence of this assumption is that the solid 
fraction must be constant over space and time. It is 
not uncommon in practical applications to violate this 
condition. This is the case namely in sedimentation, or 
when the flow of a suspension has a free surface (e.g. 
debris flow or sheet flow), or when the shear occurs at 
imposed confining stress (type II experiments). In such 
case, the divergence of the large scale velocity fields of 
both phases balance each other, giving rise to long range 
hydromechanical coupling, also known as poromechani¬ 
cal couplings in porous media theory and soil mechanics 
El . This strong coupling governs a range of phenomena 
such as liquefaction of loose materials or, conversely, de¬ 
lays in the solid-fluid transition in dense materials [25]. 
Note also that, even if both phases comove at large 
scales, the implicit fluid velocity field corresponding to 
the lubrication terms is not divergence free, and there¬ 
fore not fully consistent. 

Beside Stokesian Dynamics, numerical methods solv¬ 
ing the Stokes or Navier-Stokes equations directly for 
the interstitial fluid have been introduced [2TH36] . In all 
cases, the so-called lubrication terms are singular and 
can not be solved explicitly by computational fluid dy¬ 
namics (CFD). Capturing the divergence of these terms 
when particles approach contact would need to shrink 
the mesh to unrealistically small element sizes around 
the contact region. A practical approach is to let CFD 
compute the non-singular terms and to add the lubri¬ 
cation terms directly using analytical expressions [ 23 , 

H3|. 

In the present work, dense suspensions are simu¬ 
lated, using a particle-based method and accounting 
for three effects: contact interactions, drag forces result¬ 
ing from the poromechanical coupling, and lubrication 
forces. The contact forces and the motion of the par¬ 
ticles are computed using a Discrete Element Method 
(DEM). The poromechanical coupling is accounted for 
using the pore-scale coupling DEM-PFV (Pore Finite 
Volume) developed recently [9|8]. The DEM-PVF code 
has been extended to periodic boundary problems for 


the purpose of the present study. And finally the lu¬ 
brication forces are introduced using frame invariant 
expressions. 

One question that we will examine in this paper is 
whether accounting for the divergence free nature of 
fluid velocity field at the local scale can significantly af¬ 
fect the rheology of dense suspensions at constant vol¬ 
ume. The poromechanical coupling will be exhibited as 
a transient effect during volume changes. Another ques¬ 
tion of interest concerns the definition of the lubrica¬ 
tion forces and torques. As particles move one relative 
to one another, normal, tangential, rolling and twist¬ 
ing motions generate different effects. These effects are 
sometimes introduced selectively in numerical models, 
considering that some of them must have negligible ef¬ 
fects, a priori. Typically, only normal and shear forces 
are computed [2711241133] , and sometimes even the shear 
force is neglected l22l fll[30 ] . Hereafter, we introduce all 
possible terms in order to evaluate, a posteriori, which 
ones can be neglected. Finally, the contributions of - 
respectively - the solid contacts, the lubrication forces, 
and the poromechanical coupling will be investigated. 

2 Numerical model 

2.1 Discrete Element Model 

An explicit finite difference scheme is employed for up¬ 
dating the position of each particle in a time-marching 
algorithm. The particles move according to the New¬ 
ton’s second law. In the absence of gravitational effects, 
the motion is driven by elastic-frictional contact forces 
defined using a soft contact approach classical in the 
DEM [12] . The contact parameters are the normal and 
shear stiffnesses k n and k s , and the angle of contact 
friction qb. The contact forces are supplemented here¬ 
after with forces coming from the interstitial fluid. A 
three-dimensional implementation of the DEM as found 
in the open source software YADE is used herein. For 
more details about the implementation, please refer to 

[35]. 

2.2 DEM-PFV coupled model 

The DEM-PFV method is used to solve a pore-scale 
version of the mass balance equation which appears in 
the continuous theory of porous media and leads to the 
so-called poromechanical coupling. Only the main steps 
of the method are outlined hereafter since the details 
can be found in previous papers. Here, we assume in¬ 
compressible phases as in 9.[8] (for compressible phases 
see [29]). A tetrahedral decomposition of the pore space 
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is introduced based on regular triangulation (figure [Tj) , 
where that part of a tetrahedron occupied by the fluid 
is called a pore. From now on V denotes the volume of 
pore i. It is uniquely defined by the positions and 
sizes of the solid particles, while the rate of change V 
also depends on their velocities x^. 

An exchange of fluid between adjacent pores i and 
j is represented by the interface flux qij. The volume 
balance equation for one one pore leads to 

j= 4 

Vi - Y, (!) 

3 = 1 

Assuming a Stokes regime entails a linear relation be¬ 
tween qij and the local pressure gradient (pi — Pj)/hj 
between two pores, where kj is a reference length (see 
ED- It leads to the following relationship between pres¬ 
sure and rate of volume change: 

j =4 j =4 

Vi= ki i = T Ki i ( Pi- Pi )• ( 2 ) 

3 = 1 3=1 


In this equation Kij is the local hydraulic conductivity. 
It must reflects the small scale geometry of the pack¬ 
ing. In details, the proposed expression of depends 
on a local hydraulic radius (area of the fluid-solid 
interface divided by the fluid volume - again see El) as 


sLf&. 


Kij = a- tJ 13 


2 77 h 


( 3 ) 


where Sfj is the cross-sectional area of the pore-throat, 
77 is the viscosity of the fluid, and a can be interpreted 
as a calibration parameter, a = 1 is known to give 
good estimates of the actual permeability of glass beads 
[32] but we used a < 1 in this study. This is further 
discussed in section [3] 

Substituting V by its expression in terms of parti¬ 
cles velocity and writing equation [ 2 ] for every element 
gives a system of linear equations. In a matrix form and 
including boundary conditions, the sytem reads 


KP — Ex + Qg + Q p, 


( 4 ) 


where K is the conductivity matrix containing the terms 
P the column vector of pressure unknowns, and E 
is the matrix defining the rates of volume change of the 
elements such that Vi = (Ex)j. Q q and Q p are flux 
vectors reflecting the boundary conditions, respectively 
source terms (imposed fluxes) in Q q and imposed pres¬ 
sures in Q p . Solved at each time step, equation 0] gives 
the discrete field of fluid pressure P as function of the 
particles velocity. 

The drag forces can be deduced from the pressure 
field. They are the integrals of the pressure p and the 


viscous stress r on the surface of the particle (for a 
non-inertial fluid, the second integral can be evaluated 
based on P using momentum balance): 

Fl = / puds + / rnds (5) 

Jdr k Jdr k 

The forces Fl are linearly dependant of P, hence a 
matrix form giving the drag forces for all particles 

F / =I / P (6) 

The Fp are introduced in Newton’s second law to¬ 
gether with the forces coming from solid contacts (F c ) 
and lubrication effects (F L defined in the next section). 

I.e. 

Mx = F c + F l + I / P, (7) 

The strong two-way coupling defined by equations 
[4] and [7] is the poromechanical coupling. It is integrated 
with an explicit scheme whose accuracy has been veri¬ 
fied in [ 8 ]. 



2.3 Lubrication 

We assume that lubrication reflects the dominant dis¬ 
sipation process locally in the sheared pore fluid. The 
lubrication effects are computed between particle pairs 
as soon as they share an edge in the triangulation (fig¬ 
ure 0 ]) (the triangulation is updated dynamically during 
the simulation). They are defined for all the elementary 
motions described in figure [ 2 j We consider particles k 
and k' with radii ak and aw , linear velocities Vk and 
and angular velocities and uv, respectively. Their 
average radius is defined as a = (a^ + ) /2 and h de¬ 

notes the inter-particle distance (surface to surface). An 
arbitrary relative motion between two particles can be 
decomposed in four elementary motions corresponding 
to normal displacement (subscript n), shear displace¬ 
ment (s), rolling (r) and twisting (t). This decomposi¬ 
tion is illustrated in figure [2j In addition, we introduce 
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the angular velocity of the local frame attached to the 
interacting pair: uj n = (vk* — Vk) x n/(a^ + + h). 

Lubrication forces and torques induced by these ele¬ 
mentary motions are defined as follow: 


F „ = 


Ft = 


7TT7 


- 2(2 fl- ( 2(2 ~\~ K ) 171 


2(2 h 
h 


vt 


( 8 ) 

( 9 ) 



Normal motion Rolling motion 

Shear motion Twist motion 


r o /3 7 a 63 h . a\ r/ N , 

= tt r] a + [(a*-a,*) x n] 

( 10 ) 

C t L = 7 TT)a 2 ^ln^ [(cjfe - Wfe/) • n] n (11) 


Fig. 2 Relative motion between particles. 


or from Frankel & Acrivos US]: 

_ 7T 77 / . 7 . 7 2a + /D 

= — ( -2 a + (2a + ft) In —-— ] v t . 


where v n = ((v&/ — v&) • n) n is the normal relative 
velocity and v* = (a fe (w/e -u; n )+a fe / (w^ -a; n )) x n is an 
objective expression of the tangential relative velocity. 
In this set of equations, the normal and shear forces, F n 
and F s , are based on Frankel & Acrivos mm whereas 
C r and Ct are based on Jeffrey & Onishi flSlfTT] - the 
reason of this choice will be discussed later. The total 
lubrication force Ffc (resp. F'£,) applied by particle k’ 
on particle k (resp. by particle k on particle k') and the 
total torque Cfc (resp. C%,) applied by particle k' on 
particle k (resp. by particle k on particle k') relative to 


the particle center read: 


Fi L — — Fh — F 4- F 

r fc — r k' — r n \ r S5 

(12) 

C{ = (a k + ^)F S + C r + C t , 

(13) 

Cl = (a# + F s - C r - C t . 

(14) 


Note that F s contributes to the total torques. 

In order to check the validity of the different lu¬ 
brication approximations for different interparticle dis¬ 
tance h 3D Finite Element Method (FEM) simulations 
of Stokes flow have been carried out. Periodic bound¬ 
ary conditions were used to represent an infinite array 
of identical spheres, fixed in space but all rotating at 
the same velocity (inset of figure [3]). Figure [3] presents 
the comparison of the FEM results with that of equa¬ 
tion (9), where Ff is determined alternatively using the 
expression from Jeffrey & Onishi m- 

_ . a 

F s = nr] a In- v t , 


Both expressions are asymptotically equivalent for 
h 0. However, our results show that the second one 
is in much better agreement with the FEM result for 
small but finite distances (. h/2a < 0.1). Both expres¬ 
sions underestimate the FEM result for h/2a > 0.1, 
which is not surprising keeping in mind that they have 
been obtained from asymptotic expansions. However, 
an additional defect of the former is that it leads to 
negative torques (i.e. torques with the same sign as the 
angular velocity) for large h. This feature is unphysical 
as it leads to a net creation of energy and can severely 
alter the stability of the numerical scheme. It was thus 
concluded that the expressions of Frankel & Acrivos 
were more suitable for implementation. 

We account for the deformability of the particles 
near the contact region by combining the above nor¬ 
mal lubrication model with a linear elastic model via 
a Maxwell-type visco-elastic scheme f|23|, partly in¬ 
spired by [26]) (figure 0]) . The parameters are k n the 
contact stiffness and v n {h) is the instantaneous viscos¬ 
ity of the interaction as defined in eq. (3), such that 
= v n {h)v n . The real velocity of approach between 
the two surfaces is h = v n — , where is the elastic 

deformation given by u e n = F^/k n . Hence, the evolu¬ 
tion of the normal lubrication force obeys the differen¬ 
tial equation 


that we have to integrate over time-steps using the form 




(16) 










6 


Donia Marzougui et al. 


120 


^ c, 

O 


80 


60 


40 


20 


r\ r\ 

1- . C / \ C ; C / 

Q 

. rS.OO 


A A FEM 

— Rankel & Acrivos 

- ■ - Jeffrey & Onishi 



k 

A 

i 

\'- 

4'. 

\ \ 

A \ . 


A 

A 

A 



'A . 

\ ' 

\ \ 

\ ‘ 
s 

A . A 
























- 2 8o 


0.1 0.2 0.3 0.4 0.5 0.6 

A 

2 a 


Fig. 3 Comparaison of viscous shear torques for the case of 
rotating sphere in a regular assembly of identical particles, h 
is the surface-to-surface distance and a is the particle’s radius. 



Fig. 5 Evolution of the contact force and of the normalized 
lubrication force in the normal direction as a function of the 
gap between two particles, e defines the roughess of the par¬ 
ticles surfaces. 


3 Numerical simulations 


The same stiffness k n is used for both the DEM contact 
model and the visco-elastic lubrication model. 

Lastly, when the particles approach each other, the 
pressure in the gap tends to infinity and the normal 
lubrication force would theoretically prevent contact. 
Like many others (see e.g. |27|). we assume that contact 
will actually occur when h is of the order of the particle 
roughness 5 (figure [5]). Thus the solid contact model 
will generate a repulsive interaction even before h = 0. 
As a result of this combination between a visco-elastic 
model and roughness at contact, h will never reach 0 
practically in numerical simulations. 


Fig. 4 Visco-elastic scheme of the interaction between two 
elastic-like particles. 


3.1 Simulation setup 

The suspension is represented by a bi-periodic packing 
made of N = 1000 frictional spheres of average radius 
a = 0.025 ± 0.01 m. The physical properties are (un¬ 
less stated otherwise for sensitivity analysis) roughness 
5 = 0.035 a, density p = 2500 kg/m 3 , normal contact 
stiffness k n /a = 5 x 10 5 Pa, shear stiffness k s = k n / 2, 
and contact friction angle cp = 30°. There is no gravity. 
The numerical sample is H = 18a high, L = 12a long 
and l = 12a wide (figure [6]). It is first confined between 
two parallel plates then sheared by moving the top and 
the bottom plates at constant velocity ±V/2 = 1.5m/s. 
The boundary conditions for the top plate are the veloc¬ 
ities v x = E/2, v z = 0, the total vertical stress T y = 750 
Pa and the fluid pressure p = 0. At the bottom plate, 
v x = —E/2, v z = 0 and the fluid velocity along the y 
axis v£ = 0 (impermeable boundary). Periodic bound¬ 
ary conditions are defined along the horizontal axis for 
both the particles and the fluid. For the latest we im¬ 
pose a null pressure gradient at the macro-scale, i.e. 
VPx = VPz = 0 (see [23]). In order to avoid preferen¬ 
tial slip zones near the plates, the first layer of spheres 
in contact with a plate is fixed to the plates by highly 
cohesive contacts. We introduce the boundary stress 
vector T = F/S where F is the total force on the top 
plate and S is the horizontal cross sectional area. T y is 
constant during the deformation, while T x is a result of 
the imposed shear. 
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No ux 


Fig. 6 Simulation cell. 


3.2 Transient vs. steady state 

Figure 0 shows the evolution of the shear stress T x , the 
pressure p and the solid fraction (j) function of the defor¬ 
mation j(T) = Jq (t)dt where y(t) = V/H(t) is the 
shear rate. The numerical results are presented for two 
cases: a first case where the poromechanical coupling 
is considered and a second one without it (i.e. ignoring 
the last term in equation 0. When the shear velocity is 
applied, a transient regime is observed, characterized by 
an increase followed by a decrease of the shear stress, a 
decrease of the solid fraction and a negative pressure of 
the fluid in the coupled case. This later effect entails a 
higher effective stress in the coupled problem, explain¬ 
ing why the shear stress reaches higher values. The sys¬ 
tem evolves toward a steady state for large deforma¬ 
tions, in which the shear stress and the solid fraction 
are approximately constant and the pressure is nearly 
zero. The poromechanical coupling has no visible effect 
at steady state: the shear stress and the solid fraction 
reach similar values for both cases. 

It is to be noted that the poromechanical coupling 
entails long range effects in the system and, ultimately, 
a dependency on the problem size (H in our case). It 
is known since Terzaghi that the characteristic time of 
such process scales with pH 2 / k where k is the intrin¬ 
sic permeability. Since k scales with a 2 , the relaxation 
time of the transient regime is proportional to rj(H/a) 2 . 
A consequence is that the peak pore pressure in fig¬ 
ure 0 scales with ( H/a ) 2 , pressure gradients scale with 
H/a 2 , and finally drag forces scale with H. Important 
consequences of this feature are that the drag forces 
and the lubrication forces are not comensurable, and 
that poromechanical effects cannot be reflected as rhe¬ 
ological properties of the bulk material - it is always 
necessary to solve a coupled problem. 


Since our main focus in this study was the bulk vis¬ 
cosity of suspensions (a rheological property), we did 
not seek a fully realistic combination of drag forces and 
lubrication forces. Practically, it let us reduce the com¬ 
putation times tremendeously by setting a = 100 in 
equation [3j while it should be close to 1 for more real¬ 
istic simulations. The duration of the transient regime 
would have been multiplied by 100 approximately, and 
the timesteps of the time marching algorithm would 
have been reduced by 100, leading to an increase of the 
total simulation time by a factor 10 4 . Though qualita¬ 
tively correct, the trends seen in fig. [7] are thus quan¬ 
titatively wrong (and remember that in any case they 
are only relevant for a specific value H). 

3.3 Stress decomposition 

Besides T, a tensor representing the average stress in 
the suspension can be calculated as 

ct = <7° + ct l + pi + a 1 (17) 

cj c = y 'Ylij Fpj ® hj is the contact stress applied 
on particles in contact where hj denotes the branch 
vector between the centers of the particles i and j. 
Similarly, a L = y J2ij ^ tj ® hj is the contribution 
from lubrication forces |T, which will be further de¬ 
composed hereafter by considering separately the nor¬ 
mal and shear components of the lubrication force, p is 
the pressure associated to the poromechanical coupling 
m- v 1 = m k w k 0 Vfc reflect the inertial effects as 
defined in [2&j where mk is the mass of particle k and 
Vk is its velocity. It can be verified in figure [7| - where 
both T x and a xy are plotted - that the two expressions 
compare consistently. Hereafter, eq. (EJ> will be used 
to assess the microscale origins of the shear stress and 
their rate dependency. 

4 Results and Discussion 

4.1 Typical results 

Figure [8] includes the evolution of different terms of 
equation (H7|) (component xy) as the suspension is sheared 
at I v = 0.21. The inertial stress a^ y (not represented 
here) is negligible compared with the total stress ((J xy < 
2.5 %T X ), which indicates that the suspension is domi¬ 
nated by contacts and viscous interactions in this case 
(further discussion in the next paragraph). Second, the 
contact stress contributes to approximately half of the 
total stress (<J xy ~ 50% T x ) whereas the other half is 
due to the normal and shear lubrication forces ( a^ ~ 

30 %T X and a^ y « 20 %T X ). The different contributions 









Fig. 7 (top): The evolution of the shear stress and the solid 
fraction as a function of the deformation for two cases: with 
poromechanical coupling and without poromechanical cou¬ 
pling. (bottom): Zoom on the transient regime. 

will be further investigated for different values of the 
viscous number I v in the following. 


Fig. 8 Decomposition of the total shear stress in contact 
stress, normal lubrication stress and shear lubrication stress. 

confining pressure, fluid viscosity, and shear rate corre¬ 
sponding to a given value of I v should give the same 
result. In order to confirm this property, three series 
of simulations were conducted in which the control pa¬ 
rameters were changed independently in each series to 
produce different values of I v . The results are plotted 
in fig. [9] versus the corresponding values of I v . fi and <£> 
are nearly the same whatever the method to change I v . 
This result offers a numerical confirmation of the con¬ 
clusion of Boyer et al. j5]. This property holds only in 
the viscous regime, i.e. as long as the inertial effects can 
be neglected. As suggested in [33], this condition may 
be characterized by the ratio I/I v , where / = y/pj 2 /T y 
is the so called inertial number. Here the ratio is at most 
I/I v = 5, and it is much smaller in most cases, while 
[33] suggest I/I v ~ 10 for the transition from the vis¬ 
cous to the inertial regime. For the rest of this study it 
was decided to keep I constant and equal to 0.14, which 
leaves 77 as the only free parameter. An exception is the 
dry case (I v = 0 ) where I = 0.005. 


4.2 Viscous number 


I v is a dimensionless form of the shear rate [5], reflecting 
the magnitude of viscous effects, and is defined as: 



The key idea here is that, in the viscous regime and 
at steady state the stress ratio fi = T x /T y and the 
solid fraction </> are entirely controlled by this unique 
parameter. In other words, all possible combinations of 


4.3 Dropping terms of the hydrodynamic model 

The consequences of neglecting some of the terms de¬ 
fined by eqs. ©-CD will be examined, in terms of 
stress ratio fi and of solid fraction </>. This question is of 
interest for the developpers of numerical models, since 
many models found in the literature are not includ¬ 
ing all terms. As we have seen before, neglecting the 
poromechanical coupling (eq. [5} has a detrimental ef¬ 
fect on the transient state but has no strong effect at 
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□ changing T 


changing rj v v changing 7 



Fig. 9 Normalized shear stress and solid fraction at steady 
state versus I v . In each series the change of I v is obtained 
by changing a different parameter: confining pressure (■), 
viscosity (★), or shear rate (▼). 


isfactory stress ratio but significantly overestimate the 
solid fraction. Considering this result, one may expect 
that simulations at imposed volume (called type I in 
the introduction) and including only normal lubrica¬ 
tion forces would underestimate the shear stress even 
more than in our case. 

Further sophisticating the hydrodynamic model does 
not yield other significant changes. The rolling torques 
(yellow stars), the twist torques (green diamond) and 
the poromechanical coupling (cyan triangles) have only 
marginal effects. We recall that this conclusion holds 
at the steady state only. As discussed previously, the 
poromechanical coupling can significantly interfere in 
the transient regimes. 


4.4 Contact stress versus hydrodynamic stress 


steady state. We now focus on the different lubrication 
terms and how they affect the result at steady state. 

Figure [TO] presents the comparison of the numer¬ 
ical results with the phenomenological laws proposed 
by Boyer et al [5] on the basis of experiments. The val¬ 
ues reported in this figure have been obtained at steady 
state by increasing the fluid viscosity while keeping the 
shear rate 7 and the normal stress T y constant. Typical 
results for T x vs. 7 are shown in the inset of the fig¬ 
ure for different values of rj. Series of simulations have 
been carried out including different combinations of the 
hydrodynamic effects. Starting with the simplest case 
where only solid contacts and normal lubrication forces 
are present, the other hydrodynamic terms are added 
one by one: shear lubrication force, rolling torque, twist 
torque and drag forces due to the poromechanical cou¬ 
pling. 

The result for I v = 0 corresponds to vanishing hy¬ 
drodynamic forces. It is obviously independent of the 
hydrodynamic assumptions, and practically it is the 
result of a simulation for a dry material. The friction 
coefficient and the solid fraction obtained in this case 
closely match the values measured experimentally by 
Boyer et al. The results obtained with normal lubrica¬ 
tion forces only (blue squares) match at least qualita¬ 
tively the empirical evolution of shear stress with I v . It 
is to be noted however that the solid fraction obtained 
with this model is almost constant for I v > 0 . 02 , while 
the experiments suggest a monotonic decrease. As soon 
as the shear lubrication forces are included, the results 
(red circles) get closer to the phenomenological law for 
both /jl(I v ) and cj)(I v ). It can be concluded that shear 
lubrication forces play a key role in the rate depen¬ 
dent dilatancy, and they contribute significantly to the 
shear stress. The normal lubrication alone lead to a sat¬ 


From now on, the results that will be analyzed are ob¬ 
tained with the full model including all possible hydro- 
dynamic interactions. Figure [Tl] shows the contribution 
to the total stress of the different terms in eq. (fI7|) . The 
contact forces play a significant role for all the values 
of I v investigated. The contact stress slightly increases 
for 0 < I v < 0.1 and saturates to an almost constant 
value for larger values of I v . Lubrication stresses, both 
normal and shear components, increase almost linearly, 
and the shear stress due to the shear components is ap¬ 
proximately twice smaller than the stress coming from 
the normal components. For values of I v > 0.15 the 
sum of the two lubrication stresses exceeds the contact 
stress. This result highlights the fact that depending on 
the value of I v two regimes are observed. At low I v the 
contact interactions are dominant whereas for I v > 0.15 
the lubrication interactions dominate. This result con¬ 
firms a constitutive property inferred by Boyer et al [5]. 
Herein, the shear stress in dense suspensions is split into 
two contributions, one coming from the contacts and 
represented by the same phenomenological law as in dry 
granular flow, the other one coming from hydrodynamic 
interactions similar to a Krieger-Dougherty viscosity. 
However Boyer et al inferred this stress partition from 
macroscopic measurements. The present results allow 
to further assess the respective contributions of con¬ 
tacts and hydrodynamic interactions. Also, the choice 
of a frictional rheology for the contact stress which sat¬ 
urates for high values of I v is predicted by our discrete 
numerical simulations. We believe that this property 
is not trivial. Based on solid fraction at the largest I v 
indeed (0 ~ 0.38, far below any values that can be 
reached in dry quasi-static granular systems), one could 
expect that no solid contacts persist. Nevertheless, the 
numerical simulations confirm the proposal of Boyer et 
al. on this aspect. 
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Fig. 10 Stress ratio /r (top) and the solid fraction <j) (bottom) 
at steady state versus I v . Each symbol represents a different 
combination of lubrication terms. The solid line is the phe¬ 
nomenological law of Boyer et al. [5]. Inset: the total shear 
stress for different values of fluid viscosity. 


4.5 Roughness 

The particles roughness e appears as a key parameter in 
the lubrication model since the surface-to-surface dis¬ 
tance h may vanish as e —)> 0, a situation in wich the lu¬ 
brication terms diverge. This situation is peculiar from 
a physical point of view but it also causes major trou¬ 
bles from a numerical point of view. Arguably perfectly 



Fig. 11 The stress ratio p, and the decomposition in contact 
stress, the normal lubrication stress and the shear lubrication 
stress. The solid line is the phenomenological law of [5]. 


smooth surfaces are rare and this parameter can be jus¬ 
tified on a physical ground. It is not always clear how 
this should be accounted for in models however, and 
in our case we don’t have a precise knowledge of what 
value of e would be relevant for the spheres used by 
Boyer et al. In order to evaluate the role of this param¬ 
eter simulations were reproduced for three different val¬ 
ues of e. They are reported in figure [12J which includes 
the total stress and the contributions from contacts and 
lubrication forces. 

The lesser the roughness, the larger the contribution 
of lubrication, as one could expect. A less expected re¬ 
sult is that the contribution of contacts is not strongly 
modified. No clear trend can be distinguished, as the 
points corresponding to the smallest e may be below or 
above the others depending on the value I v . This sug¬ 
gests that the different values are simply due to impre- 
cisions in the evaluation of the steady state, and that 
the contribution of contacts could be independent of 
roughness. Overall, the difference on the total stress is 
of the order of 10% or less, showing a relatively moder¬ 
ate effect. An effect on dilatancy is visible, the smoother 
particles dilate more, but again moderate. All three val¬ 
ues seem to reflect the main features of the behaviour 
in spite of quantitative differences. Roughness does not 
appear as a key parameter here. The limit e 0 re¬ 
mains as an open question, which can’t be studied eas¬ 
ily with our method due to numerical difficulties. 
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Fig. 12 Stress ratio /r and solid fraction cf> at steady state 
for different values of the roughness parameter e. 


4.6 Microstructure 

In order to link averaged quantities to micro-scale vari¬ 
ables, we examine how various quantities depends on 
the orientation of the particle pairs (figure Based 
on the orientation of the unit normal, every interaction 
corresponds to a position on the unit sphere. For an 
arbitrary point M on the unit sphere it is possible to 
compute averages of interaction variables. For instance, 
the average distance between the spheres is defined as: 
h(M ) = J2n k edS h(k)/NM where dS is a small angular 
sector centered on M and Nm is the number of interac¬ 
tions associated to dS. In figure [13] this value is normal¬ 
ized by the particle diameter (2a). We define similarly 
the average normal velocity and average shear velocity 
normalized by 2 aj. 

The other plots of figure [13] are density functions. 
The density of contacts is obtained by counting the 
number N c of interactions with ft < e in dS , and divid¬ 
ing by the total number of spheres 7V S , so that P(M) = 
N c (M)/(dS N s ). The densities of the lubrication stress 
term ct l (M) and of the contact stress term cr c (M) are 
obtained by restricting the sums defined for eq. 
to the subset of interactions associated to M. For in- 

stance a L (M) = y } ZsT,n ij eds F ij ® 0n figure 0 
are the component xy of both stress tensors, normal¬ 
ized by the confining pressure. The lubrication stress is 
further decomposed into one part due to normal forces 
and another part due to shear forces. 

Despite the fact that all the functions introduced 
above define surfaces in the 3D space, in figure [T3l only 
the values for M in the ( Oxy ) plane are plotted. The 
results are given for three different values of I v : I v s= 0 
(i.e. a dry medium - green line in the figure), I v = 0.025 
(red line) and I v = 0.2 (blue line). All distributions are 


7r-periodic, the comments hereafter refer to the interval 
[0,7r]. 

Figure [13] a shows that there is a minimum in the 
density of contacts near 6 = 7 r /4 for all cases. Con¬ 
versely, the higher density is observed for orientations 
between 37 r /4 and 7 r. A noticeable effect of increasing 
I v is that the number of contacts near 0 = 7 r /4 vanish. 
This peculiar effect makes a clear difference between 
the viscous regime and the dry regime as P{0 = 7 r/ 4 ) 
is always strictly positive in the latest. 

The average distance (figure [J3]b) is increasing with 
I v on overall. Since increasing I v corresponds to a de¬ 
crease of solid fraction, this trend is not surprising. The 
average distance is anisotropic and takes larger values 
near 6 = 7 r/ 4 , consistently with the lower density of 
contacts in this region. 

The normal component of normalized relative ve¬ 
locity (figure [13] c) is positive on [0, 7r/2] and negative 
on [7r/2,7r]. The extrema are of the same order in both 
cases, although it can be noted that the velocity of ap¬ 
proaching particles ( [7t/ 2, 7r] ) is slightly lower. The shear 
component (figure [13] d) is positive on - approximately 
- [7r/4,37r/4] and negative on [0, 7r/4] U [37t/4, tt\. In this 
case the extrema are clearly different, with the largest 
relative velocities near 0 = it/ 2. We note that the graph 
of the shear velocity is not exactly symmetric as the 
maximum values are in fact a bit before 6 = 7r/2. It 
can be explained by smaller ft and more solid contacts 
preventing sliding when 0 > 7r/2. 

The magnitude of the normalized relative velocity is 
increasing slightly with I v . This is observed for both the 
normal and the shear components. This trend can be 
explained by considering the growing average distance 
between particles. If all particles were simply following 
the mean field velocity, then the relative velocity would 
obviously grow with ft. Since the local fluctuations of 
velocity with respect to the mean field are not modi¬ 
fying the relative velocities in average, this correlation 
holds. 

The lubrication forces for a given relative velocity 
are decreasing functions of ft. Thus, it could be antici¬ 
pated from figure [T3]b that the contribution of viscous 
interactions to the bulk stress is dominated by interac¬ 
tions oriented along 0 = 37 r/ 4 , which have in average 
smaller values of ft but nearly similar values of normal 
velocity. The results of figure [13] e show a quite different 
picture. The density of stress due to the normal lubri¬ 
cation forces is actually slightly larger on [ 0 , ir/2] and it 
matches the shape of the normal velocity closely (with 
the difference that the sign is always positive due to 
the branch l in the diadic product ®Z). This feature 
may be explained by the strong correlations between ft 
and v n . 
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The contribution of the shear lubrication forces to 
the bulk stress is dominated by the interactions near 
6 = 7r/2, consistently with the evolution of average 
shear velocity (the fact that the contribution vanishes 
when 9 = 0 is an effect of the product 0 l with l 
nearly horizontal). Like v s , the density of stress reaches 
a pick slightly before 0 = tt/2. 

Considering the cumulated contributions cr L = cr^-\- 
o-g, the non-symmetry of the density of stress becomes 
even more visible. The larger contribution is due to in¬ 
teractions in the range [7r/4, n/2]. This is also the range 
in which there are nearly no solid contacts between the 
particles. Thus the non-symmetry may be due to a com¬ 
plex interplay between the solid contacts and the vis¬ 
cous interactions. 

The density of contact stress (figure [I3jf) is strongly 
anisotropic, as expected from the lack of contacts on 
[0,7r/2]. When I v = 0 some contacts in [0,7r/2] con¬ 
tribute negatively, but the contribution is small. For 
I v > 0 this contribution becomes negligible. On over¬ 
all, the contacts contribute to the bulk stress mainly via 
repulsive forces in the direction near 0 = 37r/4. The con¬ 
tribution of contacts to the bulk stress increases with 
I v . This feature was already observed in figure [10] on 
the macroscopic variables. 


5 Conclusion 

We described a complete modeling framework for the 
numerical simulation of dense suspensions. The model 
includes solid contacts between particles using a soft 
contact approach, short range hydrodynamic interac¬ 
tions defined by frame-invariant expressions of forces 
and torques in the lubrication approximation, and the 
poromechanical coupling solved using the DEM-PFV 
technique. 

Numerical experiments of simple shear at imposed 
confining normal stress have been conducted in an at¬ 
tempt to reproduce recent rheometer experiments on 
beads. The simulations are in excellent agreement with 
the empirical data in terms of bulk shear stress and 
solid fraction at the steady state, for a range of dimen¬ 
sionless shear rate corresponding to 0 < I v <0.45. 

The poromechanical coupling was shown to have a 
significant effect in the transient regime when the de¬ 
formation starts. However, no significant effects of this 
coupling have been exhibited at steady state. 

The results obtained by neglecting some of the lubri¬ 
cation terms leads to the following conclusions. First, 
the normal lubrication term has the larger contribu¬ 
tion to the bulk stress. However, considering this term 
alone leads to underestimate the shear stress, and with 


such simplification the model is unable to reflect the 
change of solid fraction for increasing I v . Combining 
both normal and shear lubrication terms gives much 
better results in terms of stress and solid fraction. Fur¬ 
ther sophistication of the model by including the terms 
associated with rolling and twisting gives only marginal 
improvements. 

The analysis of the various contributions to the bulk 
stress: contact forces, hydrodynamic forces and fluid 
pressure, has lead to the following conclusions. The con¬ 
tribution of contacts to the bulk shear stress in the 
permanent regime increases with I v . This result may 
be seen as counter-intuitive. First, as higher I v leads to 
lower solid fraction, one would expect contacts contri¬ 
bution to decrease progressively with I v and ultimately 
vanish. Second, the assumption that the lubrication ef¬ 
fects would prevent solid contacts in suspensions has 
been the cornerstone of many theoretical and numerical 
models in the past. Instead, the simulations reported in 
the present paper suggest that both the contact stress 
and the lubrication stress increase monotonically in the 
range of I v investigated. 

As roughness of particules is decreased, the contri¬ 
bution of lubrication forces is increased to some extent. 
However, it does not lead to a decrease of contact forces, 
which remain nearly unchanged. Again, it is against the 
idea that lubrication is preventing solid contacts. The 
reasoning leading to this idea may be simply flawn due 
to conceptual mistakes. Namely, lubrication forces are 
often perceived as repulsive force wereas they are neu¬ 
tral overall, inhibiting the opening and the closure of 
contacts almost equally as revealed by microstructural 
variables. 

The distribution of micro-structural variables re¬ 
vealed a complex interplay between the contact fabric 
and the hydrodynamic interactions. The anisotropy of 
contact orientation appears to be more pronounced in 
suspensions as compared to dry granular materials, due 
to the effect of the hydrodynamic interactions. This can 
explain at least partly why the contribution of contacts 
to the bulk shear stress is increasing. 

This contact stress may reach a maximum and ul¬ 
timately vanish for larger I v . When does it happen re¬ 
mains an open question. Capturing this transition in 
numerical simulations is a great challenge for numerical 
models since the lubrication terms do not reflect appro¬ 
priately all hydrodynamic interactions in more dilute 
regimes. 
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Fig. 13 Distributions of normalized quantities in the (x,y) plane [top/left] for the dry case, I v = 0.025 and I v = 0.2: PDF 
of contact orientation, normal velocity, shear velocity, lubrication stress which is the sum of the normal and shear lubrication 
stress, and contact stress. 
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